Overview

Bayesian distribution regression (BDR) is built on a few pieces of theory:

  1. Kernel mean embeddings
  2. Gaussian process regression (as the prior for the kernel mean embeddings)
  3. Landmark points for dimensionality reduction (advantages of this over FastFood? both are reducing the dimensionality of the feature embedding so that we can calualate the empirical mean)

This notebook builds the theory of Bayesian Distribution Regression piece by piece:

  1. Explore kernels from the kernlab package. We’ll later use the functions from this package to calculate mean embeddings for bags of observations \(\{x_j^i\}_{j=1}^{N_i}\).

Problem statement

We observe survey responses \((\mathbf{x}_i, y_i)\) for \(i \in S\), where \(\mathbf{x}_i \in \mathbb{R}^{V_s}\) is a vector of \(v \in (1, \dots V_s)\) covariates (i.e. age, gender, education, income) and \(y_i\) denotes who a respondent says they’re likely to support in an upcoming political election. \(S \subset (1, \dots, N)\) denotes the set of indicies in the population of size \(N\) that were observed in the survey.

We also observe a much larger set of \(V_u >> V_s\) covariates \(\mathbf{x}_i \in \mathbb{R}^{V_u}\), \(V_s \subset V_u\), for all \(i \in D\), where \(D \subset N\). The covariates \(v \in V_u \setminus V_s\) are highly predictive of \(y_i\).

We say that a survey respondent is “matched” to the population if \(i \in M\) where \(M = S \cap D\).

The goal here is to estimate \(\hat{y}_i \forall i \in D\). Usually to do this, predictive models are built using only the “matched” data, or using respondents \(i \in M\). However, we would like to leverage \(i \in S \setminus M\), for which we observe \(y_i\) and \(x_{i, 1}, \dots, x_{i, V_s}\) but not \(\mathbf{x}_{i, \mathbf{v}}\) for \(\mathbf{v} \in V_u \setminus V_s\).

Question:

Explore Kernels from kernlab package

## create a RBF kernel function with sigma hyper-parameter 0.05
rbf = rbfdot(sigma = 1)
## create artificial data set
x <- matrix(rnorm(60), 6, 10)
## compute kernel matrix
kx <- kernelMatrix(rbf, x)  ## k_12 is equivalent to exp(- sum((x[1,] - x[2,])^2))

Basic Gaussian process

A Gaussian process is a collection of random variables, any finite subset of which follows a multivariate normal distribution

prior: \(f \sim \mathcal{GP}(\mathbf{0}, k(x, x'))\)

We fix a finite set of \(d\) points \(\mathbf{s} = (s_1, \dots, s_d)\) in order to draw from the GP as a multivariate normal. Conditioning on \(s\), this becomes \(\mathbf{f} \sim \mathcal{N}(\mathbf{0}_d, \Sigma_d)\) where \(\Sigma_d\) is the empirical covariance matrix of our points \(s\).

We then ``observe" 5 points \(\{x_i, y_i\}_{i = 1}^n\). We assume that y is observed without noise such that \(y = f(x)\). In order to predict at a new point \(x_\star\), we condiditon on the observed points \(X\) and their associated outcomes \(Y\), and arrive at the closed-form posterior predictive distribution of \(f_\star\)

posterior: \(f_\star | x_\star, X, \mathbf{y} \sim \mathcal{N}(k_\star K^{-1}\mathbf{y}, k_{\star\star} - k_\star K^{-1}k_\star^T)\)

where \(K_{ij} = k(x_i, x_j)\) and \(k_\star = [k(x_1, x_\star), \dots, k(x_n, x_\star)]\) and \(k_{\star\star} = k(x_star, x_star)\)

# define function to calculate empirical covariance matrix based on Gaussian kernel with length-scale l
calcSigma = function(x1, x2, l = 1){
  Sigma = matrix(nrow = length(x1), ncol = length(x2))
  
  for(i in 1:length(x1)){
    for(j in 1:length(x2)){
      Sigma[i,j] = exp(-1/2 * ((x1[i] - x2[j])/l)^2)
    }
  }
  return(Sigma)
}
# define the points that we want to fix for the function draw
s = seq(-5, 5, by = 0.1)
# calculate the empirical covariance of the points (basically 0)
Sigma = calcSigma(s, s)
# draw 5 samples from the GP prior at the fixed points
x_prior = t(mvrnorm(n = 5, mu = rep(0, length(s)), Sigma = Sigma))
x_draws = data.frame(cbind(s, x_prior))
x_draws = melt(x_draws, id = 's', value.name = 'f(x)')
# plot draws from the prior
plot_gp_prior = ggplot(x_draws, aes(x = s, y = `f(x)`, color = variable)) + geom_line() +
  ggtitle("Draws from GP Prior") + theme_minimal()
# specify points that we observe
x_obs = c(-4, -2.5, -1, 0, 4)
y_obs = c(-2, 0, 1, 2, -1)
# calculate kernel matricies needed for posterior evaluation
K_xx = calcSigma(x_obs, x_obs)
K_xstarx = calcSigma(s, x_obs)
K_xxstar = calcSigma(x_obs, s)
K_xstarxstar = calcSigma(s, s)
# calculate posterior mean and variance
mu = K_xstarx %*% solve(K_xx) %*% y_obs
Sigma_star = K_xstarxstar - K_xstarx %*% solve(K_xx) %*% K_xxstar
# draw 5 samples from the posterior at fixed eval points s
x_post = t(mvrnorm(5, mu = mu, Sigma = Sigma_star))
x_post = data.frame(cbind(s, x_post))
x_post = melt(x_post, id = 's', value.name = 'f_star')
# calculate 95% confidence bands for posterior
f_star_covar = data.frame(cbind(s, ub = mu + 2 * sqrt(diag(Sigma_star)), lb = mu - 2 * sqrt(diag(Sigma_star))))
NaNs producedNaNs produced
# plot draws from posterior
plot_gp_post = ggplot() + 
  geom_ribbon(data = f_star_covar, aes(x = s, ymin = V3, ymax = V2), fill = "grey", alpha = 0.5) +
  geom_line(data = x_post, aes(x = s, y = f_star, color = variable)) +
  #geom_point(aes(x = x_obs, y = y_obs)) +
  ggtitle("Draws from GP Posterior") +
  theme_minimal()
  
grid.arrange(plot_gp_prior, plot_gp_post, nrow = 1)

Try out BDR on survey data

Goal: Estimate individual-level support using BDR on state-level support outcomes

Data

Undeclared level(s) 2, 3, 4 added in variable: densityUndeclared level(s) 2, 3, 4 added in variable: sdensityUndeclared level(s) 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95 added in variable: ageUndeclared level(s) 1, 2, 3, 4, 5, 6, 7 added in variable: hh1Undeclared level(s) 1, 2, 3, 4, 5, 6, 7 added in variable: hh3

Individual-level survey responses (n = nrow(data_sept18) from September 2018 Pew Research survey (https://www.people-press.org/dataset/september-2018-political-survey/).

Recode survey responses

Recode the survey responses - combine q7 (strong) and q8 (leaners) to get full support response

## support 
data_sept18[, .N, .(q7, q8)]
data_sept18[, qsupport := NULL]
Adding new column 'qsupport' then assigning NULL (deleting it).
data_sept18[q7 == "Democratic Party's candidate" | q8 == "Democratic Party's candidate", qsupport := '1-D']
data_sept18[q7 == "Republican Party's candidate" | q8 == "Republican Party's candidate", qsupport := '2-R']
data_sept18[q8 == "(VOL) Other", qsupport := '3-O']
data_sept18[is.na(qsupport), qsupport := '4-DK/R']
data_sept18[, .(.N, pct = .N/nrow(data_sept18)), .(q7, q8, qsupport)][order(qsupport)]

calculate percentages by state

data_state = data_sept18[, .(y_dem = mean(qsupport == '1-D')
                             , y_rep = mean(qsupport == '2-R')
                             , y_other = mean(qsupport == '3-O')
                             , y_dkr = mean(qsupport == '4-DK/R')
                             , total_respondents = .N
                             # sample covariates
                             , age_under30 = mean(as.numeric(age) < 30)
                             , race_W = mean(racecmb == 'White')
                             ), by = sstate][order(sstate)]
head(data_state)
data_sept18[, .N, hh1][order(hh1)]
X = data_sept18[, .(state = as.character(sstate)
                    , sex_male = as.numeric(sex == 'Male')
                    , sex_female = as.numeric(sex == 'Female')
                    , age = as.numeric(age)
                    
                    , educ_postgrad = as.numeric(grepl("Postgraduate", educ))
                    , educ_bach = as.numeric(grepl("Four year", educ))
                    , educ_assoc = as.numeric(grepl("Two year associate|Sone college", educ))
                    , educ_highschool = as.numeric(grepl("High school graduate", educ))
                    , educ_none = as.numeric(grepl("Less than high school|High school incomplete", educ))
                    
                    , hisp = as.numeric(hisp == 'Yes')
                    , race_white = as.numeric(racecmb == 'White')
                    , race_black = as.numeric(racecmb == 'Black or African-American')
                    , race_asian = as.numeric(racecmb == 'Asian or Asian-American')
                    , race_mixedother = as.numeric(racecmb == 'Mixed Race' | racecmb == 'Or some other race')
                    
                    , relig_protestant = as.numeric(grepl("Protestant|Christian", relig))
                    , relig_catholic = as.numeric(grepl("Catholic", relig))
                    , relig_athiest = as.numeric(grepl("Athiest|Agnostic", relig))
                    , relig_jewish = as.numeric(grepl("Jewish", relig))
                    , relig_muslim = as.numeric(grepl("Muslim", relig))
                    , relig_LDS = as.numeric(grepl("Mormon", relig))
                    
                    , hh_n = ifelse(hh1 == '8 or more', 8, ifelse(hh1 == "Don't know/Refused", 2, as.numeric(hh1)))
                    
                    , income_under10 = as.numeric(income == 'Less than $10,000')
                    , income_10to20 = as.numeric(income == '10 to under $20,000')
                    , income_20to30 = as.numeric(income == '20 to under $30,000')
                    , income_30to40 = as.numeric(income == '30 to under $40,000')
                    , income_40to50 = as.numeric(income == '40 to under $50,000')
                    , income_50to75 = as.numeric(income == '50 to under $75,000')
                    , income_75to100 = as.numeric(income == '75 to under $100,000')
                    , income_100to150 = as.numeric(income == '100 to under $150,000')
                    , income_over150 = as.numeric(income == 'Over $150,000')
                    , income_refused = as.numeric(income == "(VOL) Don't know/Refused")
                    )]
# center and scale X
X_numeric = scale(X[, 2:ncol(X)])
# replace NA and NaN with 0
X_numeric = apply(X_numeric, 2, function(x) {
  x[is.na(x) | is.nan(x)] <- 0
  x
})
# create full covar matrix
X = data.frame(state = X$state, X_numeric)
##### Create test/training split
set.seed(123)
# train_ind = createDataPartition(X
#                                 , y = X$state
#                                 , p = .6, list = FALSE)
train_ind = sample.int(nrow(X), size = 1000, replace = F)
X_train = X[train_ind, ]
X_test = X[-train_ind, ]
Y_train = data_sept18[train_ind, .(state = as.character(sstate), y_dem = as.numeric(qsupport == '1-D'))]
Y_test = data_sept18[-train_ind, .(state = as.character(sstate), y_dem = as.numeric(qsupport == '1-D'))]
# create data that we actually observe
Y_train_agg = Y_train[, .(y_dem_pct = mean(y_dem), .N), by = state]
Y_test_agg = Y_test[, .(y_dem_pct = mean(y_dem), .N), by = state]
Y_train_agg

Mean embeddings

We observe bags of points \(\{x_j^i\}_{j = 1}^{N_i}\) from \(i = 1, \dots, n\), where \(x_j^i \in \mathbb{R}^p\). We have to embed the \(x_j^i\) in feature space and then take the mean. This means we can’t use the kernel trick which would take us directly from \(\mathbb{R}^p \times \mathbb{R}^p \to \mathbb{R}\). In order to calculate the coordinate of each observation in feature space, we convert each point to an explicit featurization: \(x_j^i \in \mathbb{R}\) is mapped to \[\phi(x_j^i) = [k(x_j^i, u_1), \dots, k(x_j^i, u_d)]^T \in \mathbb{R}^d\] where \(\mathbf{u} = \{u_l\}_{l = 1}^d\) are a set of landmark points. For now we’ll set \(d = 100\) and choose \(u_l\) using k-means clustering.

Specify landmark points \(\mathbf{u} = (u_1, \dots, u_k)\) as the centroids of k-means clustering

n_centers = 100
groups = kmeans(X_train[, -1], centers = n_centers)
u = as.matrix(groups$centers)
head(u)
     sex_male  sex_female          age educ_postgrad   educ_bach educ_assoc educ_highschool  educ_none       hisp race_white
1  0.90084435 -0.90084435  0.006395155   -0.35356593  1.65478397 -0.3815408      -0.4857439 -0.2324988 -0.3766247  0.3196829
2  0.23075063 -0.23075063  0.302619660   -0.35356593 -0.60396396 -0.3815408      -0.4857439  4.2986449 -0.3766247  0.5618264
3  0.09673189 -0.09673189 -0.314968403   -0.03553791  0.07366042  0.2186588      -0.4857439 -0.2324988  1.4415397 -1.7788945
4  0.90084435 -0.90084435 -0.718954580    2.82671430 -0.60396396 -0.3815408      -0.4857439 -0.2324988 -0.3766247  0.5618264
5  0.90084435 -0.90084435 -0.806020567   -0.35356593  0.52541001 -0.3815408      -0.4857439  2.0330730  1.1385123 -1.7788945
6 -1.10943680  1.10943680  0.587035216   -0.35356593 -0.15221437 -0.3815408       0.7011146 -0.2324988 -0.3766247  0.5618264
  race_black race_asian race_mixedother relig_protestant relig_catholic relig_athiest relig_jewish relig_muslim  relig_LDS
1 -0.2660408 -0.1546641      -0.3474161      -0.89153704      0.8273477    -0.2112912   -0.1658855  -0.08297413 -0.1447158
2 -0.2660408 -0.1546641      -0.3474161       0.67378425     -0.4866399    -0.2112912   -0.1658855  -0.08297413 -0.1447158
3 -0.2660408 -0.1546641       2.8767521      -0.08651466      0.7835481    -0.2112912   -0.1658855  -0.08297413 -0.1447158
4 -0.2660408 -0.1546641      -0.3474161       1.12101890     -0.4866399    -0.2112912   -0.1658855  -0.08297413 -0.1447158
5 -0.2660408  3.1536379       1.2646680      -0.89153704     -0.4866399    -0.2112912   -0.1658855  -0.08297413  6.9061577
6 -0.2660408 -0.1546641      -0.3474161       0.58433732     -0.1479231    -0.2112912   -0.1658855  -0.08297413 -0.1447158
         hh_n income_under10 income_10to20 income_20to30 income_30to40 income_40to50 income_50to75 income_75to100
1  0.20112327     -0.2053351    -0.2745375    -0.3101534    -0.2967138     -0.261095    -0.3815408    -0.36469037
2  0.02109454     -0.2053351    -0.2745375     0.8673554     1.3321941     -0.261095    -0.3815408    -0.01967046
3  0.04379381     -0.2053351    -0.2745375    -0.3101534    -0.2967138     -0.261095    -0.3815408     2.74048880
4  1.38305102     -0.2053351    -0.2745375    -0.3101534    -0.2967138     -0.261095    -0.3815408     2.74048880
5  1.95053289     -0.2053351    -0.2745375     1.4561097     1.5358076     -0.261095    -0.3815408    -0.36469037
6 -0.56908661     -0.2053351    -0.2745375     3.2223729    -0.2967138     -0.261095    -0.3815408    -0.36469037
  income_100to150 income_over150 income_refused
1               0              0     -0.2600654
2               0              0     -0.3666934
3               0              0     -0.3666934
4               0              0     -0.3666934
5               0              0     -0.3666934
6               0              0     -0.3666934

Create RBF kernel and get embeddings of training set \(\phi(x)\) using kernel and landmark points \[\phi(x_i) = [k(x_i, u_1), \dots, k(x_i, u_k)]\]

Choose scale parameter \(\sigma\) for RBF kernel - use median heuristic for now

rbf1 = rbfdot(sigma = 1)
K_sigma = kernelMatrix(rbf1, x = as.matrix(X_train[, -1]), y = u)
sigma = median(K_sigma)
ggplot(melt(K_sigma), aes(x = log(value))) + geom_histogram() + geom_vline(xintercept = log(sigma), color = 'red')

rbf = rbfdot(sigma = 0.02)
phi_x = kernelMatrix(rbf, x = as.matrix(X_train[, -1]), y = u)
phi_x = data.table(state = X_train$state, phi_x)
setnames(phi_x, c('state', paste0('u', 1:n_centers)))
head(phi_x)

Calculate mean embeddings for each state

mu_hat = phi_x[, lapply(.SD, mean), .SDcols = names(phi_x)[2:ncol(phi_x)], by = state][order(state)]
mu_hat

Model 1 - Frequentist

Normal LASSO regression on the kernel mean embeddings - no Bayesian treatment.

\[y_i = \beta^T\hat{\mu}_i + b\] Doesn’t account for uncertainty in the coefficients \(\beta\) or in the mean embeddings \(\hat{\mu}_i\), or observation error

Build the model

Use 5-fold CV kernel LASSO to estimate \(\beta\), \(b\) for \(\hat{y}_i = \beta^T\hat{\mu}_i + b\). Use the first LASSO to select non-zero coefs, then re-fit without penalty to accont for covariate shrinkage.

# fit initial lasso to find significant covars
fit_lambda = cv.glmnet(x = as.matrix(mu_hat[, -1, with = F])
                       , y = Y_train_agg[order(state), ]$y_dem_pct
                       , nfolds = 7)

# get indicies of non-0 covars (except intercept)
ind = which(coef(fit_lambda, s = 'lambda.min')[-1] != 0)

# re-fit LASSO only with non-0 covars
fit = glmnet(x = as.matrix(mu_hat[, ind + 1, with = F])
             , y = Y_train_agg[order(state), ]$y_dem_pct)

Predict

Get embeddings of test points and empirical means of bags of embedded test points. Use same landmark points \(u\) as for the training set.

# get embeddings of test points
phi_x_test = kernelMatrix(rbf, x = as.matrix(X_test[, -1]), y = u)
phi_x_test = data.table(state = X_test$state, phi_x_test)
setnames(phi_x_test, c('state', paste0('u', 1:n_centers)))

# calculate empirical means
mu_hat_test = phi_x_test[, lapply(.SD, mean), .SDcols = names(phi_x_test)[2:ncol(phi_x_test)], by = state][order(state)]
mu_hat_test_gender = cbind(phi_x_test, sex_female = as.numeric(X_test$sex_female > 0))[, lapply(.SD, mean), .SDcols = names(phi_x_test)[2:ncol(phi_x_test)], by = .(state, sex_female)][order(state, sex_female)]

First, test predictions at the state-level. We can also test at the state-gender level.

y_hat_agg = predict(fit, newx = as.matrix(mu_hat_test[, ind + 1, with = F]), s = 0)
y_hat_agg = cbind(Y_test[, .(y_dem_pct = mean(y_dem), .N), by = state][order(state)], y_hat_agg)

ggplot(y_hat_agg, aes(x = y_dem_pct, y = `1`, size = N)) + 
  geom_point() + 
  ylab("Predicted % Dem") +
  xlab("Actual % Dem") +
  ggtitle("Predicted by actual % Dem support across states")

y_hat_agg_gender = predict(fit, newx = as.matrix(mu_hat_test_gender[, ind + 1, with = F]), s = 0)
y_hat_agg_gender = cbind(cbind(Y_test, sex_female = as.numeric(X_test$sex_female > 0))[, .(y_dem_pct = mean(y_dem), .N), by = .(state, sex_female)][order(state, sex_female)], y_hat_agg_gender)

ggplot(y_hat_agg_gender, aes(x = y_dem_pct, y = `1`, size = N, color = as.factor(sex_female))) + 
  geom_point() + 
  geom_abline(intercept = 0, slope = 1) +
  ylab("Predicted % Dem") +
  xlab("Actual % Dem") +
  ggtitle("Predicted by actual % Dem support across states & genders")

print(paste0("MSE (state): ", mean((y_hat_agg[,`1`] - y_hat_agg$y_dem_pct)^2)))
print(paste0("MSE (state & gender): ", mean((y_hat_agg_gender[,`1`] - y_hat_agg_gender$y_dem_pct)^2)))

Predict at the individual-level

# predict on new (individual-level) data set of X_test embedded in feature space
y_hat = predict(fit, newx = as.matrix(phi_x_test[, ind + 1, with = F]), s= 0)  #no regularization
y_hat = cbind(Y_test, y_hat)

# calculate overall classification rate
class_rates = y_hat[, .(.N
          , pct_correct = mean((`1` > 0.5 & y_dem == 1) | (`1` <= 0.5 & y_dem == 0))
          , pct_false_neg = mean(`1` <= 0.5 & y_dem != 0)
          , pct_false_pos = mean(`1` > 0.5 & y_dem != 1)
          , mse = mean((y_dem - `1`)^2)
          )]
class_rates

Plot actual percent Dem by decile of the score

# get deciles of predicted score
y_hat[, pred_decile := cut(`1`, breaks = quantile(`1`, probs = seq(0,1,0.1)), labels = 1:10, include.lowest = T)]

#ggplot(y_hat, aes(x = `1`, color = as.factor(y_dem)), alpha = 0.2) + geom_density()

#plot actual pct of each decile that supports dems
ggplot(y_hat[, .(pct_y_dem = mean(y_dem)), by = pred_decile], aes(x = pred_decile, y = pct_y_dem)) + 
  geom_bar(stat = 'identity') +
  ggtitle("Pct Dem supporter by score decile")

Model 2 - Bayesian Linear Regression

Adding uncertainty: incorporate uncertainty of the coefficients \(\beta\)

We can interpret this in 2 ways: 1) Bayesian linear regression 2) Gaussian process regression conditioned on finite set of observed points.

  1. Bayesian linear regression \[\beta \sim \mathcal{N}(0, \Sigma_p)\\ y_i | \beta, \mathbf{x}_i \sim \mathcal{N}(\beta^T\hat{\mu}_i, \sigma^2)\]

For disitribution regression, we substitute in the empirical mean embeddings \(\hat{\mu}_i\) for \(x_i\).

  1. GP regression conditioned on observed points (weight space view) \[ f(\mathbf{x}) = \phi(\mathbf{x})^T\mathbf{w}\\ y = f(\mathbf{x}) + \epsilon\\ \epsilon \sim \mathcal{N}(0, \sigma_n^2)\\ \mathbf{w} \sim \mathcal{N}(0, \Sigma_p) \]

we can equivalently write \[ f \sim \mathcal{GP}(0, k(x, x'))\\ y_i|f(x_i) \sim \mathcal{N}(0, K + I\sigma^2) \] where \(k(x,x')\) is the kernel corresponding to the feature map \(\phi\)

For GPR, the feature map \(\phi(x_i)\) corresponds to the kernel that was used to embed \(x_i\) in feature space. Using landmark points \(\mathbf{u} = \{u_l\}_{l=1}^d\), this becomes \(\phi(x_i) = [k(x_i, u_1), \dots, k(x_i, u_l)]\)

We can interpret the weights \(\mathbf{w}\) of GP regression as the coefficients \(\beta\) in BLR. \(\Sigma_p\) is the prior covariance matrix for the weights/coefficients.

Fit GP

Now, we predict the GP at the bag-level


y_hat_agg = predict(fit, newx = as.matrix(mu_hat_test[, ind + 1, with = F]), s = 0)
y_hat_agg = cbind(Y_test[, .(y_dem_pct = mean(y_dem), .N), by = state][order(state)], y_hat_agg)

ggplot(y_hat_agg, aes(x = y_dem_pct, y = `1`, size = N)) + 
  geom_point() + 
  ylab("Predicted % Dem") +
  xlab("Actual % Dem") +
  ggtitle("Predicted by actual % Dem support across states")

And at the individual level

y_hat_gp = predict(fit_gp, newdata = phi_x_test[, -1])
y_hat_gp = cbind(Y_test, y_hat_gp)

y_hat_gp[, pred_decile := cut(V1, breaks = quantile(V1, probs = seq(0,1,0.1)), labels = 1:10, include.lowest = T)]

#plot actual pct of each decile that supports dems
ggplot(y_hat_gp[, .(pct_y_dem = mean(y_dem)), by = pred_decile], aes(x = pred_decile, y = pct_y_dem)) + 
  geom_bar(stat = 'identity') +
  ggtitle("Pct Dem supporter by score decile")

THINGS TO DO

“Unlabeled data can simply be used as centers” - Szabo 2016 - not entirely applicable here, we observe labels, just at the aggregate level

General questions

  • how much to focus on this specific example v. the more general question of “does dist reg work for individuals”?
  • What are the reasons that you’re skeptical that BDR will work for individuals? What are the shortcomings that you envision?

  • advantages to implementing base model as GPR v. BLR? GPR more flexible for non-Gaussian outcomes (i.e. multinomial)? Can combine kernels (i.e. mean embeddings + spatial component)? Downside is computation?

  • batching?

Choose data set:

  • census microdata with actual county-level election outcomes. pro: large (like actual voterfile), con: outcomes not at the individual-level
  • Pew survey data. pro: outcomes at the individual-level, con: it’s small
  • census microdata for covar data with survey data for outomes. Pro: large covar data, closer to real scenario, con: can only link at the state-level, still no individual-level outcome data
  • synthetic data. pro: can make it exactly what we want, con: effort to synthesize

** benefit of the larger data set - more realistic, we’d really need the landmark points when we don’t if only working with the survey data

** larger data also has more covariate data available

parameter tuning:

  • number of landmark points (how is this normally chosen?) - Leon’s code has 30 as the default for 1000 test, 1000 training (similar scale to data here)
  • Calculate landmark points with full data set or just training?
  • scale parameter \(\sigma\) for RBF kernel (right now just usiing median heuristic)

Choosing bag definitions

  • Since we have individual-level labels in the survey data, we can choose how we aggregate it (do we even need to aggregate it? YES because we need the strata definitions in order to link it back to the full file)
  • Bags based on demos instead of state (more specific to this poll than most we’d be dealing with)?

RBF Network v. Single RBF kernel

  • relationship between neural network implementation (re Leon’s code) and the kernel mean embedding? Just that we can exactly interpret kernel embedding as a neural network and maybe easier to implement/fast optimization of marginal log likelihood?

Benchmark:

  • Compare against lasso/custom-built logit or maybe NN that has access to the individual-level outcomes (more realistic benchmark than), just using a biased subset of the training observations (mimic asymmetric matching)
  • MRP - standard for EI in political applications
  • straight stratification using covariates observed in the survey

Implement full Bayesian treatment (uncertainty in mean embeddings as well as bag size)

  • flesh out full generative model

Structured observational uncertainty

  • heriarchical model?
LS0tCnRpdGxlOiAiQmF5ZXNpYW4gREIgUHJhY3RpY2UiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCgpgYGB7ciBzZXR1cCwgZWNobz1GfQpybShsaXN0ID0gbHMoKSkKbGlicmFyeShkYXRhLnRhYmxlKQpsaWJyYXJ5KGZvcmVpZ24pCmxpYnJhcnkoa2VybmxhYikKbGlicmFyeShNQVNTKSAgIyBmb3IgbXZybm9ybQpsaWJyYXJ5KGdncGxvdDIpCmxpYnJhcnkoZ3JpZEV4dHJhKQpsaWJyYXJ5KGdsbW5ldCkKbGlicmFyeShjYXJldCkKc2V0d2QoJ34vZ2l0aHViL2Jkci8nKQpgYGAKCiMjIE92ZXJ2aWV3CgpCYXllc2lhbiBkaXN0cmlidXRpb24gcmVncmVzc2lvbiAoQkRSKSBpcyBidWlsdCBvbiBhIGZldyBwaWVjZXMgb2YgdGhlb3J5OgoKMS4gS2VybmVsIG1lYW4gZW1iZWRkaW5ncwoyLiBHYXVzc2lhbiBwcm9jZXNzIHJlZ3Jlc3Npb24gKGFzIHRoZSBwcmlvciBmb3IgdGhlIGtlcm5lbCBtZWFuIGVtYmVkZGluZ3MpCjMuIExhbmRtYXJrIHBvaW50cyBmb3IgZGltZW5zaW9uYWxpdHkgcmVkdWN0aW9uIChhZHZhbnRhZ2VzIG9mIHRoaXMgb3ZlciBGYXN0Rm9vZD8gYm90aCBhcmUgcmVkdWNpbmcgdGhlIGRpbWVuc2lvbmFsaXR5IG9mIHRoZSBmZWF0dXJlIGVtYmVkZGluZyBzbyB0aGF0IHdlIGNhbiBjYWx1YWxhdGUgdGhlIGVtcGlyaWNhbCBtZWFuKQoKVGhpcyBub3RlYm9vayBidWlsZHMgdGhlIHRoZW9yeSBvZiBCYXllc2lhbiBEaXN0cmlidXRpb24gUmVncmVzc2lvbiBwaWVjZSBieSBwaWVjZToKCjEuIEV4cGxvcmUga2VybmVscyBmcm9tIHRoZSBga2VybmxhYmAgcGFja2FnZS4gIFdlJ2xsIGxhdGVyIHVzZSB0aGUgZnVuY3Rpb25zIGZyb20gdGhpcyBwYWNrYWdlIHRvIGNhbGN1bGF0ZSBtZWFuIGVtYmVkZGluZ3MgZm9yIGJhZ3Mgb2Ygb2JzZXJ2YXRpb25zICRce3hfal5pXH1fe2o9MX1ee05faX0kLgoKIyMgUHJvYmxlbSBzdGF0ZW1lbnQKCldlIG9ic2VydmUgc3VydmV5IHJlc3BvbnNlcyAkKFxtYXRoYmZ7eH1faSwgeV9pKSQgZm9yICRpIFxpbiBTJCwgd2hlcmUgJFxtYXRoYmZ7eH1faSBcaW4gXG1hdGhiYntSfV57Vl9zfSQgaXMgYSB2ZWN0b3Igb2YgJHYgXGluICgxLCBcZG90cyBWX3MpJCBjb3ZhcmlhdGVzIChpLmUuIGFnZSwgZ2VuZGVyLCBlZHVjYXRpb24sIGluY29tZSkgYW5kICR5X2kkIGRlbm90ZXMgd2hvIGEgcmVzcG9uZGVudCBzYXlzIHRoZXkncmUgbGlrZWx5IHRvIHN1cHBvcnQgaW4gYW4gdXBjb21pbmcgcG9saXRpY2FsIGVsZWN0aW9uLiAkUyBcc3Vic2V0ICgxLCBcZG90cywgTikkIGRlbm90ZXMgdGhlIHNldCBvZiBpbmRpY2llcyBpbiB0aGUgcG9wdWxhdGlvbiBvZiBzaXplICROJCB0aGF0IHdlcmUgb2JzZXJ2ZWQgaW4gdGhlIHN1cnZleS4gIAoKV2UgYWxzbyBvYnNlcnZlIGEgbXVjaCBsYXJnZXIgc2V0IG9mICRWX3UgPj4gVl9zJCBjb3ZhcmlhdGVzICRcbWF0aGJme3h9X2kgXGluIFxtYXRoYmJ7Un1ee1ZfdX0kLCAkVl9zIFxzdWJzZXQgVl91JCwgZm9yIGFsbCAkaSBcaW4gRCQsIHdoZXJlICREIFxzdWJzZXQgTiQuICBUaGUgY292YXJpYXRlcyAkdiBcaW4gVl91IFxzZXRtaW51cyBWX3MkIGFyZSBoaWdobHkgcHJlZGljdGl2ZSBvZiAkeV9pJC4KCldlIHNheSB0aGF0IGEgc3VydmV5IHJlc3BvbmRlbnQgaXMgIm1hdGNoZWQiIHRvIHRoZSBwb3B1bGF0aW9uIGlmICRpIFxpbiBNJCB3aGVyZSAkTSA9IFMgXGNhcCBEJC4gIAoKVGhlIGdvYWwgaGVyZSBpcyB0byBlc3RpbWF0ZSAkXGhhdHt5fV9pIFxmb3JhbGwgaSBcaW4gRCQuICBVc3VhbGx5IHRvIGRvIHRoaXMsIHByZWRpY3RpdmUgbW9kZWxzIGFyZSBidWlsdCB1c2luZyBvbmx5IHRoZSAibWF0Y2hlZCIgZGF0YSwgb3IgdXNpbmcgcmVzcG9uZGVudHMgJGkgXGluIE0kLiAgSG93ZXZlciwgd2Ugd291bGQgbGlrZSB0byBsZXZlcmFnZSAkaSBcaW4gUyBcc2V0bWludXMgTSQsIGZvciB3aGljaCB3ZSBvYnNlcnZlICR5X2kkIGFuZCAkeF97aSwgMX0sIFxkb3RzLCB4X3tpLCBWX3N9JCBidXQgbm90ICRcbWF0aGJme3h9X3tpLCBcbWF0aGJme3Z9fSQgZm9yICRcbWF0aGJme3Z9IFxpbiBWX3UgXHNldG1pbnVzIFZfcyQuCgoKClF1ZXN0aW9uOgoKLSBJcyB0aGlzIGJldHRlciB0aGFuIHByZWRpY3RpbmcgdXNpbmcgdGhlIGNsYXNzIHByb3BvcnRpb24gaW4gZWFjaCBzdHJhdHVtIGRlZmluZWQgYnkgdGhlIGNvdmFyaWF0ZXMgb2JzZXJ2ZWQgaW4gdGhlIHNhbXBsZSAkcF9zJD8KCgoKIyMgRXhwbG9yZSBLZXJuZWxzIGZyb20gYGtlcm5sYWJgIHBhY2thZ2UKCmBgYHtyIGtlcm5lbHN9CiMjIGNyZWF0ZSBhIFJCRiBrZXJuZWwgZnVuY3Rpb24gd2l0aCBzaWdtYSBoeXBlci1wYXJhbWV0ZXIgMC4wNQpyYmYgPSByYmZkb3Qoc2lnbWEgPSAxKQoKIyMgY3JlYXRlIGFydGlmaWNpYWwgZGF0YSBzZXQKeCA8LSBtYXRyaXgocm5vcm0oNjApLCA2LCAxMCkKCiMjIGNvbXB1dGUga2VybmVsIG1hdHJpeApreCA8LSBrZXJuZWxNYXRyaXgocmJmLCB4KSAgIyMga18xMiBpcyBlcXVpdmFsZW50IHRvIGV4cCgtIHN1bSgoeFsxLF0gLSB4WzIsXSleMikpCmBgYAoKCiMjIEJhc2ljIEdhdXNzaWFuIHByb2Nlc3MKCkEgR2F1c3NpYW4gcHJvY2VzcyBpcyBhIGNvbGxlY3Rpb24gb2YgcmFuZG9tIHZhcmlhYmxlcywgYW55IGZpbml0ZSBzdWJzZXQgb2Ygd2hpY2ggZm9sbG93cyBhIG11bHRpdmFyaWF0ZSBub3JtYWwgZGlzdHJpYnV0aW9uCgoqcHJpb3IqOiAkZiBcc2ltIFxtYXRoY2Fse0dQfShcbWF0aGJmezB9LCBrKHgsIHgnKSkkCgpXZSBmaXggYSBmaW5pdGUgc2V0IG9mICRkJCBwb2ludHMgJFxtYXRoYmZ7c30gPSAoc18xLCBcZG90cywgc19kKSQgaW4gb3JkZXIgdG8gZHJhdyBmcm9tIHRoZSBHUCBhcyBhIG11bHRpdmFyaWF0ZSBub3JtYWwuIENvbmRpdGlvbmluZyBvbiAkcyQsIHRoaXMgYmVjb21lcyAkXG1hdGhiZntmfSBcc2ltIFxtYXRoY2Fse059KFxtYXRoYmZ7MH1fZCwgXFNpZ21hX2QpJCB3aGVyZSAkXFNpZ21hX2QkIGlzIHRoZSBlbXBpcmljYWwgY292YXJpYW5jZSBtYXRyaXggb2Ygb3VyIHBvaW50cyAkcyQuCgpXZSB0aGVuIGBgb2JzZXJ2ZSIgNSBwb2ludHMgJFx7eF9pLCB5X2lcfV97aSA9IDF9Xm4kLiAgV2UgYXNzdW1lIHRoYXQgeSBpcyBvYnNlcnZlZCB3aXRob3V0IG5vaXNlIHN1Y2ggdGhhdCAkeSA9IGYoeCkkLiAgSW4gb3JkZXIgdG8gcHJlZGljdCBhdCBhIG5ldyBwb2ludCAkeF9cc3RhciQsIHdlIGNvbmRpZGl0b24gb24gdGhlIG9ic2VydmVkIHBvaW50cyAkWCQgYW5kIHRoZWlyIGFzc29jaWF0ZWQgb3V0Y29tZXMgJFkkLCBhbmQgYXJyaXZlIGF0IHRoZSBjbG9zZWQtZm9ybSBwb3N0ZXJpb3IgcHJlZGljdGl2ZSBkaXN0cmlidXRpb24gb2YgJGZfXHN0YXIkCgoqcG9zdGVyaW9yKjogJGZfXHN0YXIgfCB4X1xzdGFyLCBYLCBcbWF0aGJme3l9IFxzaW0gXG1hdGhjYWx7Tn0oa19cc3RhciBLXnstMX1cbWF0aGJme3l9LCBrX3tcc3RhclxzdGFyfSAtIGtfXHN0YXIgS157LTF9a19cc3Rhcl5UKSQKCndoZXJlICRLX3tpan0gPSBrKHhfaSwgeF9qKSQgYW5kICRrX1xzdGFyID0gW2soeF8xLCB4X1xzdGFyKSwgXGRvdHMsIGsoeF9uLCB4X1xzdGFyKV0kIGFuZCAka197XHN0YXJcc3Rhcn0gPSBrKHhfc3RhciwgeF9zdGFyKSQKCmBgYHtyIHNpbXBsZS1HUCwgZmlnLmhlaWdodD0zLCBmaWcud2lkdGg9Nn0KIyBkZWZpbmUgZnVuY3Rpb24gdG8gY2FsY3VsYXRlIGVtcGlyaWNhbCBjb3ZhcmlhbmNlIG1hdHJpeCBiYXNlZCBvbiBHYXVzc2lhbiBrZXJuZWwgd2l0aCBsZW5ndGgtc2NhbGUgbApjYWxjU2lnbWEgPSBmdW5jdGlvbih4MSwgeDIsIGwgPSAxKXsKICBTaWdtYSA9IG1hdHJpeChucm93ID0gbGVuZ3RoKHgxKSwgbmNvbCA9IGxlbmd0aCh4MikpCiAgCiAgZm9yKGkgaW4gMTpsZW5ndGgoeDEpKXsKICAgIGZvcihqIGluIDE6bGVuZ3RoKHgyKSl7CiAgICAgIFNpZ21hW2ksal0gPSBleHAoLTEvMiAqICgoeDFbaV0gLSB4MltqXSkvbCleMikKICAgIH0KICB9CiAgcmV0dXJuKFNpZ21hKQp9CgojIGRlZmluZSB0aGUgcG9pbnRzIHRoYXQgd2Ugd2FudCB0byBmaXggZm9yIHRoZSBmdW5jdGlvbiBkcmF3CnMgPSBzZXEoLTUsIDUsIGJ5ID0gMC4xKQoKIyBjYWxjdWxhdGUgdGhlIGVtcGlyaWNhbCBjb3ZhcmlhbmNlIG9mIHRoZSBwb2ludHMgKGJhc2ljYWxseSAwKQpTaWdtYSA9IGNhbGNTaWdtYShzLCBzKQoKIyBkcmF3IDUgc2FtcGxlcyBmcm9tIHRoZSBHUCBwcmlvciBhdCB0aGUgZml4ZWQgcG9pbnRzCnhfcHJpb3IgPSB0KG12cm5vcm0obiA9IDUsIG11ID0gcmVwKDAsIGxlbmd0aChzKSksIFNpZ21hID0gU2lnbWEpKQp4X2RyYXdzID0gZGF0YS5mcmFtZShjYmluZChzLCB4X3ByaW9yKSkKeF9kcmF3cyA9IG1lbHQoeF9kcmF3cywgaWQgPSAncycsIHZhbHVlLm5hbWUgPSAnZih4KScpCgojIHBsb3QgZHJhd3MgZnJvbSB0aGUgcHJpb3IKcGxvdF9ncF9wcmlvciA9IGdncGxvdCh4X2RyYXdzLCBhZXMoeCA9IHMsIHkgPSBgZih4KWAsIGNvbG9yID0gdmFyaWFibGUpKSArIGdlb21fbGluZSgpICsKICBnZ3RpdGxlKCJEcmF3cyBmcm9tIEdQIFByaW9yIikgKyB0aGVtZV9taW5pbWFsKCkKCiMgc3BlY2lmeSBwb2ludHMgdGhhdCB3ZSBvYnNlcnZlCnhfb2JzID0gYygtNCwgLTIuNSwgLTEsIDAsIDQpCnlfb2JzID0gYygtMiwgMCwgMSwgMiwgLTEpCgojIGNhbGN1bGF0ZSBrZXJuZWwgbWF0cmljaWVzIG5lZWRlZCBmb3IgcG9zdGVyaW9yIGV2YWx1YXRpb24KS194eCA9IGNhbGNTaWdtYSh4X29icywgeF9vYnMpCktfeHN0YXJ4ID0gY2FsY1NpZ21hKHMsIHhfb2JzKQpLX3h4c3RhciA9IGNhbGNTaWdtYSh4X29icywgcykKS194c3RhcnhzdGFyID0gY2FsY1NpZ21hKHMsIHMpCgojIGNhbGN1bGF0ZSBwb3N0ZXJpb3IgbWVhbiBhbmQgdmFyaWFuY2UKbXUgPSBLX3hzdGFyeCAlKiUgc29sdmUoS194eCkgJSolIHlfb2JzClNpZ21hX3N0YXIgPSBLX3hzdGFyeHN0YXIgLSBLX3hzdGFyeCAlKiUgc29sdmUoS194eCkgJSolIEtfeHhzdGFyCgojIGRyYXcgNSBzYW1wbGVzIGZyb20gdGhlIHBvc3RlcmlvciBhdCBmaXhlZCBldmFsIHBvaW50cyBzCnhfcG9zdCA9IHQobXZybm9ybSg1LCBtdSA9IG11LCBTaWdtYSA9IFNpZ21hX3N0YXIpKQp4X3Bvc3QgPSBkYXRhLmZyYW1lKGNiaW5kKHMsIHhfcG9zdCkpCnhfcG9zdCA9IG1lbHQoeF9wb3N0LCBpZCA9ICdzJywgdmFsdWUubmFtZSA9ICdmX3N0YXInKQoKIyBjYWxjdWxhdGUgOTUlIGNvbmZpZGVuY2UgYmFuZHMgZm9yIHBvc3RlcmlvcgpmX3N0YXJfY292YXIgPSBkYXRhLmZyYW1lKGNiaW5kKHMsIHViID0gbXUgKyAyICogc3FydChkaWFnKFNpZ21hX3N0YXIpKSwgbGIgPSBtdSAtIDIgKiBzcXJ0KGRpYWcoU2lnbWFfc3RhcikpKSkKCiMgcGxvdCBkcmF3cyBmcm9tIHBvc3RlcmlvcgpwbG90X2dwX3Bvc3QgPSBnZ3Bsb3QoKSArIAogIGdlb21fcmliYm9uKGRhdGEgPSBmX3N0YXJfY292YXIsIGFlcyh4ID0gcywgeW1pbiA9IFYzLCB5bWF4ID0gVjIpLCBmaWxsID0gImdyZXkiLCBhbHBoYSA9IDAuNSkgKwogIGdlb21fbGluZShkYXRhID0geF9wb3N0LCBhZXMoeCA9IHMsIHkgPSBmX3N0YXIsIGNvbG9yID0gdmFyaWFibGUpKSArCiAgI2dlb21fcG9pbnQoYWVzKHggPSB4X29icywgeSA9IHlfb2JzKSkgKwogIGdndGl0bGUoIkRyYXdzIGZyb20gR1AgUG9zdGVyaW9yIikgKwogIHRoZW1lX21pbmltYWwoKQogIApncmlkLmFycmFuZ2UocGxvdF9ncF9wcmlvciwgcGxvdF9ncF9wb3N0LCBucm93ID0gMSkKYGBgCgoKIyMgVHJ5IG91dCBCRFIgb24gc3VydmV5IGRhdGEKCipHb2FsKjogRXN0aW1hdGUgaW5kaXZpZHVhbC1sZXZlbCBzdXBwb3J0IHVzaW5nIEJEUiBvbiBzdGF0ZS1sZXZlbCBzdXBwb3J0IG91dGNvbWVzIAoKIyMjIERhdGEKCmBgYHtyIGxvYWQtZGF0YSwgZWNobyA9IEYsIG1lc3NhZ2U9IEZ9CgojIG5vdGU6IGdlbyB2YXJpYWJsZXMgdGhhdCBiZWdpbiB3aXRoICdzJyBhcmUgc2VsZi1yZXBvcnRlZCAoSSB0aGluayBiYXNlZCBvbiB0aGUgcmVhZG1lKSAtIG1vcmUgYWNjdXJhdGUKCmRhdGFfc2VwdDE4ID0gZGF0YS50YWJsZShyZWFkLnNwc3MoJ2RhdGEvU2VwdDE4L1NlcHQxOCBwdWJsaWMuc2F2JywgdG8uZGF0YS5mcmFtZSA9IFQpLCBzdHJpbmdzQXNGYWN0b3JzID0gRikKZGF0YV9zZXB0MTgKYGBgCgpJbmRpdmlkdWFsLWxldmVsIHN1cnZleSByZXNwb25zZXMgKG4gPSBgbnJvdyhkYXRhX3NlcHQxOGApIGZyb20gU2VwdGVtYmVyIDIwMTggUGV3IFJlc2VhcmNoIHN1cnZleSAoaHR0cHM6Ly93d3cucGVvcGxlLXByZXNzLm9yZy9kYXRhc2V0L3NlcHRlbWJlci0yMDE4LXBvbGl0aWNhbC1zdXJ2ZXkvKS4KCgoKIyMjIyBSZWNvZGUgc3VydmV5IHJlc3BvbnNlcwoKUmVjb2RlIHRoZSBzdXJ2ZXkgcmVzcG9uc2VzIC0gY29tYmluZSBxNyAoc3Ryb25nKSBhbmQgcTggKGxlYW5lcnMpIHRvIGdldCBmdWxsIHN1cHBvcnQgcmVzcG9uc2UKYGBge3IgcmVjb2RlLWRhdGEsIG1lc3NhZ2UgPSBGfQojIyBzdXBwb3J0IApkYXRhX3NlcHQxOFssIC5OLCAuKHE3LCBxOCldCmRhdGFfc2VwdDE4WywgcXN1cHBvcnQgOj0gTlVMTF0KZGF0YV9zZXB0MThbcTcgPT0gIkRlbW9jcmF0aWMgUGFydHkncyBjYW5kaWRhdGUiIHwgcTggPT0gIkRlbW9jcmF0aWMgUGFydHkncyBjYW5kaWRhdGUiLCBxc3VwcG9ydCA6PSAnMS1EJ10KZGF0YV9zZXB0MThbcTcgPT0gIlJlcHVibGljYW4gUGFydHkncyBjYW5kaWRhdGUiIHwgcTggPT0gIlJlcHVibGljYW4gUGFydHkncyBjYW5kaWRhdGUiLCBxc3VwcG9ydCA6PSAnMi1SJ10KZGF0YV9zZXB0MThbcTggPT0gIihWT0wpIE90aGVyIiwgcXN1cHBvcnQgOj0gJzMtTyddCmRhdGFfc2VwdDE4W2lzLm5hKHFzdXBwb3J0KSwgcXN1cHBvcnQgOj0gJzQtREsvUiddCgpkYXRhX3NlcHQxOFssIC4oLk4sIHBjdCA9IC5OL25yb3coZGF0YV9zZXB0MTgpKSwgLihxNywgcTgsIHFzdXBwb3J0KV1bb3JkZXIocXN1cHBvcnQpXQpgYGAKCiMgY2FsY3VsYXRlIHBlcmNlbnRhZ2VzIGJ5IHN0YXRlCmBgYHtyfQpkYXRhX3N0YXRlID0gZGF0YV9zZXB0MThbLCAuKHlfZGVtID0gbWVhbihxc3VwcG9ydCA9PSAnMS1EJykKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAsIHlfcmVwID0gbWVhbihxc3VwcG9ydCA9PSAnMi1SJykKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAsIHlfb3RoZXIgPSBtZWFuKHFzdXBwb3J0ID09ICczLU8nKQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICwgeV9ka3IgPSBtZWFuKHFzdXBwb3J0ID09ICc0LURLL1InKQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICwgdG90YWxfcmVzcG9uZGVudHMgPSAuTgogICAgICAgICAgICAgICAgICAgICAgICAgICAgICMgc2FtcGxlIGNvdmFyaWF0ZXMKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAsIGFnZV91bmRlcjMwID0gbWVhbihhcy5udW1lcmljKGFnZSkgPCAzMCkKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAsIHJhY2VfVyA9IG1lYW4ocmFjZWNtYiA9PSAnV2hpdGUnKQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICksIGJ5ID0gc3N0YXRlXVtvcmRlcihzc3RhdGUpXQoKaGVhZChkYXRhX3N0YXRlKQpgYGAKCgpgYGB7cn0KZGF0YV9zZXB0MThbLCAuTiwgaGgxXVtvcmRlcihoaDEpXQoKWCA9IGRhdGFfc2VwdDE4WywgLihzdGF0ZSA9IGFzLmNoYXJhY3Rlcihzc3RhdGUpCiAgICAgICAgICAgICAgICAgICAgLCBzZXhfbWFsZSA9IGFzLm51bWVyaWMoc2V4ID09ICdNYWxlJykKICAgICAgICAgICAgICAgICAgICAsIHNleF9mZW1hbGUgPSBhcy5udW1lcmljKHNleCA9PSAnRmVtYWxlJykKICAgICAgICAgICAgICAgICAgICAsIGFnZSA9IGFzLm51bWVyaWMoYWdlKQogICAgICAgICAgICAgICAgICAgIAogICAgICAgICAgICAgICAgICAgICwgZWR1Y19wb3N0Z3JhZCA9IGFzLm51bWVyaWMoZ3JlcGwoIlBvc3RncmFkdWF0ZSIsIGVkdWMpKQogICAgICAgICAgICAgICAgICAgICwgZWR1Y19iYWNoID0gYXMubnVtZXJpYyhncmVwbCgiRm91ciB5ZWFyIiwgZWR1YykpCiAgICAgICAgICAgICAgICAgICAgLCBlZHVjX2Fzc29jID0gYXMubnVtZXJpYyhncmVwbCgiVHdvIHllYXIgYXNzb2NpYXRlfFNvbmUgY29sbGVnZSIsIGVkdWMpKQogICAgICAgICAgICAgICAgICAgICwgZWR1Y19oaWdoc2Nob29sID0gYXMubnVtZXJpYyhncmVwbCgiSGlnaCBzY2hvb2wgZ3JhZHVhdGUiLCBlZHVjKSkKICAgICAgICAgICAgICAgICAgICAsIGVkdWNfbm9uZSA9IGFzLm51bWVyaWMoZ3JlcGwoIkxlc3MgdGhhbiBoaWdoIHNjaG9vbHxIaWdoIHNjaG9vbCBpbmNvbXBsZXRlIiwgZWR1YykpCiAgICAgICAgICAgICAgICAgICAgCiAgICAgICAgICAgICAgICAgICAgLCBoaXNwID0gYXMubnVtZXJpYyhoaXNwID09ICdZZXMnKQogICAgICAgICAgICAgICAgICAgICwgcmFjZV93aGl0ZSA9IGFzLm51bWVyaWMocmFjZWNtYiA9PSAnV2hpdGUnKQogICAgICAgICAgICAgICAgICAgICwgcmFjZV9ibGFjayA9IGFzLm51bWVyaWMocmFjZWNtYiA9PSAnQmxhY2sgb3IgQWZyaWNhbi1BbWVyaWNhbicpCiAgICAgICAgICAgICAgICAgICAgLCByYWNlX2FzaWFuID0gYXMubnVtZXJpYyhyYWNlY21iID09ICdBc2lhbiBvciBBc2lhbi1BbWVyaWNhbicpCiAgICAgICAgICAgICAgICAgICAgLCByYWNlX21peGVkb3RoZXIgPSBhcy5udW1lcmljKHJhY2VjbWIgPT0gJ01peGVkIFJhY2UnIHwgcmFjZWNtYiA9PSAnT3Igc29tZSBvdGhlciByYWNlJykKICAgICAgICAgICAgICAgICAgICAKICAgICAgICAgICAgICAgICAgICAsIHJlbGlnX3Byb3Rlc3RhbnQgPSBhcy5udW1lcmljKGdyZXBsKCJQcm90ZXN0YW50fENocmlzdGlhbiIsIHJlbGlnKSkKICAgICAgICAgICAgICAgICAgICAsIHJlbGlnX2NhdGhvbGljID0gYXMubnVtZXJpYyhncmVwbCgiQ2F0aG9saWMiLCByZWxpZykpCiAgICAgICAgICAgICAgICAgICAgLCByZWxpZ19hdGhpZXN0ID0gYXMubnVtZXJpYyhncmVwbCgiQXRoaWVzdHxBZ25vc3RpYyIsIHJlbGlnKSkKICAgICAgICAgICAgICAgICAgICAsIHJlbGlnX2pld2lzaCA9IGFzLm51bWVyaWMoZ3JlcGwoIkpld2lzaCIsIHJlbGlnKSkKICAgICAgICAgICAgICAgICAgICAsIHJlbGlnX211c2xpbSA9IGFzLm51bWVyaWMoZ3JlcGwoIk11c2xpbSIsIHJlbGlnKSkKICAgICAgICAgICAgICAgICAgICAsIHJlbGlnX0xEUyA9IGFzLm51bWVyaWMoZ3JlcGwoIk1vcm1vbiIsIHJlbGlnKSkKICAgICAgICAgICAgICAgICAgICAKICAgICAgICAgICAgICAgICAgICAsIGhoX24gPSBpZmVsc2UoaGgxID09ICc4IG9yIG1vcmUnLCA4LCBpZmVsc2UoaGgxID09ICJEb24ndCBrbm93L1JlZnVzZWQiLCAyLCBhcy5udW1lcmljKGhoMSkpKQogICAgICAgICAgICAgICAgICAgIAogICAgICAgICAgICAgICAgICAgICwgaW5jb21lX3VuZGVyMTAgPSBhcy5udW1lcmljKGluY29tZSA9PSAnTGVzcyB0aGFuICQxMCwwMDAnKQogICAgICAgICAgICAgICAgICAgICwgaW5jb21lXzEwdG8yMCA9IGFzLm51bWVyaWMoaW5jb21lID09ICcxMCB0byB1bmRlciAkMjAsMDAwJykKICAgICAgICAgICAgICAgICAgICAsIGluY29tZV8yMHRvMzAgPSBhcy5udW1lcmljKGluY29tZSA9PSAnMjAgdG8gdW5kZXIgJDMwLDAwMCcpCiAgICAgICAgICAgICAgICAgICAgLCBpbmNvbWVfMzB0bzQwID0gYXMubnVtZXJpYyhpbmNvbWUgPT0gJzMwIHRvIHVuZGVyICQ0MCwwMDAnKQogICAgICAgICAgICAgICAgICAgICwgaW5jb21lXzQwdG81MCA9IGFzLm51bWVyaWMoaW5jb21lID09ICc0MCB0byB1bmRlciAkNTAsMDAwJykKICAgICAgICAgICAgICAgICAgICAsIGluY29tZV81MHRvNzUgPSBhcy5udW1lcmljKGluY29tZSA9PSAnNTAgdG8gdW5kZXIgJDc1LDAwMCcpCiAgICAgICAgICAgICAgICAgICAgLCBpbmNvbWVfNzV0bzEwMCA9IGFzLm51bWVyaWMoaW5jb21lID09ICc3NSB0byB1bmRlciAkMTAwLDAwMCcpCiAgICAgICAgICAgICAgICAgICAgLCBpbmNvbWVfMTAwdG8xNTAgPSBhcy5udW1lcmljKGluY29tZSA9PSAnMTAwIHRvIHVuZGVyICQxNTAsMDAwJykKICAgICAgICAgICAgICAgICAgICAsIGluY29tZV9vdmVyMTUwID0gYXMubnVtZXJpYyhpbmNvbWUgPT0gJ092ZXIgJDE1MCwwMDAnKQogICAgICAgICAgICAgICAgICAgICwgaW5jb21lX3JlZnVzZWQgPSBhcy5udW1lcmljKGluY29tZSA9PSAiKFZPTCkgRG9uJ3Qga25vdy9SZWZ1c2VkIikKICAgICAgICAgICAgICAgICAgICApXQoKIyBjZW50ZXIgYW5kIHNjYWxlIFgKWF9udW1lcmljID0gc2NhbGUoWFssIDI6bmNvbChYKV0pCgojIHJlcGxhY2UgTkEgYW5kIE5hTiB3aXRoIDAKWF9udW1lcmljID0gYXBwbHkoWF9udW1lcmljLCAyLCBmdW5jdGlvbih4KSB7CiAgeFtpcy5uYSh4KSB8IGlzLm5hbih4KV0gPC0gMAogIHgKfSkKIyBjcmVhdGUgZnVsbCBjb3ZhciBtYXRyaXgKWCA9IGRhdGEuZnJhbWUoc3RhdGUgPSBYJHN0YXRlLCBYX251bWVyaWMpCgoKIyMjIyMgQ3JlYXRlIHRlc3QvdHJhaW5pbmcgc3BsaXQKc2V0LnNlZWQoMTIzKQojIHRyYWluX2luZCA9IGNyZWF0ZURhdGFQYXJ0aXRpb24oWAojICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgLCB5ID0gWCRzdGF0ZQojICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgLCBwID0gLjYsIGxpc3QgPSBGQUxTRSkKCnRyYWluX2luZCA9IHNhbXBsZS5pbnQobnJvdyhYKSwgc2l6ZSA9IDEwMDAsIHJlcGxhY2UgPSBGKQoKWF90cmFpbiA9IFhbdHJhaW5faW5kLCBdClhfdGVzdCA9IFhbLXRyYWluX2luZCwgXQoKWV90cmFpbiA9IGRhdGFfc2VwdDE4W3RyYWluX2luZCwgLihzdGF0ZSA9IGFzLmNoYXJhY3Rlcihzc3RhdGUpLCB5X2RlbSA9IGFzLm51bWVyaWMocXN1cHBvcnQgPT0gJzEtRCcpKV0KWV90ZXN0ID0gZGF0YV9zZXB0MThbLXRyYWluX2luZCwgLihzdGF0ZSA9IGFzLmNoYXJhY3Rlcihzc3RhdGUpLCB5X2RlbSA9IGFzLm51bWVyaWMocXN1cHBvcnQgPT0gJzEtRCcpKV0KCiMgY3JlYXRlIGRhdGEgdGhhdCB3ZSBhY3R1YWxseSBvYnNlcnZlCllfdHJhaW5fYWdnID0gWV90cmFpblssIC4oeV9kZW1fcGN0ID0gbWVhbih5X2RlbSksIC5OKSwgYnkgPSBzdGF0ZV0KWV90ZXN0X2FnZyA9IFlfdGVzdFssIC4oeV9kZW1fcGN0ID0gbWVhbih5X2RlbSksIC5OKSwgYnkgPSBzdGF0ZV0KCllfdHJhaW5fYWdnCmBgYAoKIyMgTWVhbiBlbWJlZGRpbmdzCgpXZSBvYnNlcnZlIGJhZ3Mgb2YgcG9pbnRzICRce3hfal5pXH1fe2ogPSAxfV57Tl9pfSQgZnJvbSAkaSA9IDEsIFxkb3RzLCBuJCwgd2hlcmUgJHhfal5pIFxpbiBcbWF0aGJie1J9XnAkLiAgV2UgaGF2ZSB0byBlbWJlZCB0aGUgJHhfal5pJCBpbiBmZWF0dXJlIHNwYWNlIGFuZCB0aGVuIHRha2UgdGhlIG1lYW4uICBUaGlzIG1lYW5zIHdlIGNhbid0IHVzZSB0aGUga2VybmVsIHRyaWNrIHdoaWNoIHdvdWxkIHRha2UgdXMgZGlyZWN0bHkgZnJvbSAkXG1hdGhiYntSfV5wIFx0aW1lcyBcbWF0aGJie1J9XnAgXHRvIFxtYXRoYmJ7Un0kLiAgSW4gb3JkZXIgdG8gY2FsY3VsYXRlIHRoZSBjb29yZGluYXRlIG9mIGVhY2ggb2JzZXJ2YXRpb24gaW4gZmVhdHVyZSBzcGFjZSwgd2UgY29udmVydCBlYWNoIHBvaW50IHRvIGFuICpleHBsaWNpdCogZmVhdHVyaXphdGlvbjogICR4X2peaSBcaW4gXG1hdGhiYntSfSQgaXMgbWFwcGVkIHRvICQkXHBoaSh4X2peaSkgPSBbayh4X2peaSwgdV8xKSwgXGRvdHMsIGsoeF9qXmksIHVfZCldXlQgXGluIFxtYXRoYmJ7Un1eZCQkCndoZXJlICRcbWF0aGJme3V9ID0gXHt1X2xcfV97bCA9IDF9XmQkIGFyZSBhIHNldCBvZiBsYW5kbWFyayBwb2ludHMuICBGb3Igbm93IHdlJ2xsIHNldCAkZCA9IDEwMCQgYW5kIGNob29zZSAkdV9sJCB1c2luZyBrLW1lYW5zIGNsdXN0ZXJpbmcuCgoKU3BlY2lmeSBsYW5kbWFyayBwb2ludHMgJFxtYXRoYmZ7dX0gPSAodV8xLCBcZG90cywgdV9rKSQgYXMgdGhlIGNlbnRyb2lkcyBvZiBrLW1lYW5zIGNsdXN0ZXJpbmcKCmBgYHtyfQpuX2NlbnRlcnMgPSAxMDAKZ3JvdXBzID0ga21lYW5zKFhfdHJhaW5bLCAtMV0sIGNlbnRlcnMgPSBuX2NlbnRlcnMpCnUgPSBhcy5tYXRyaXgoZ3JvdXBzJGNlbnRlcnMpCmhlYWQodSkKYGBgCgoKQ3JlYXRlIFJCRiBrZXJuZWwgYW5kIGdldCBlbWJlZGRpbmdzIG9mIHRyYWluaW5nIHNldCAkXHBoaSh4KSQgdXNpbmcga2VybmVsIGFuZCBsYW5kbWFyayBwb2ludHMKJCRccGhpKHhfaSkgPSBbayh4X2ksIHVfMSksIFxkb3RzLCBrKHhfaSwgdV9rKV0kJAoKQ2hvb3NlIHNjYWxlIHBhcmFtZXRlciAkXHNpZ21hJCBmb3IgUkJGIGtlcm5lbCAtIHVzZSBtZWRpYW4gaGV1cmlzdGljIGZvciBub3cKYGBge3J9CnJiZjEgPSByYmZkb3Qoc2lnbWEgPSAxKQoKS19zaWdtYSA9IGtlcm5lbE1hdHJpeChyYmYxLCB4ID0gYXMubWF0cml4KFhfdHJhaW5bLCAtMV0pLCB5ID0gdSkKc2lnbWEgPSBtZWRpYW4oS19zaWdtYSkKCmdncGxvdChtZWx0KEtfc2lnbWEpLCBhZXMoeCA9IGxvZyh2YWx1ZSkpKSArIGdlb21faGlzdG9ncmFtKCkgKyBnZW9tX3ZsaW5lKHhpbnRlcmNlcHQgPSBsb2coc2lnbWEpLCBjb2xvciA9ICdyZWQnKQpgYGAKCgpgYGB7cn0KcmJmID0gcmJmZG90KHNpZ21hID0gMC4wMikKcGhpX3ggPSBrZXJuZWxNYXRyaXgocmJmLCB4ID0gYXMubWF0cml4KFhfdHJhaW5bLCAtMV0pLCB5ID0gdSkKcGhpX3ggPSBkYXRhLnRhYmxlKHN0YXRlID0gWF90cmFpbiRzdGF0ZSwgcGhpX3gpCnNldG5hbWVzKHBoaV94LCBjKCdzdGF0ZScsIHBhc3RlMCgndScsIDE6bl9jZW50ZXJzKSkpCgpoZWFkKHBoaV94KQpgYGAKCkNhbGN1bGF0ZSBtZWFuIGVtYmVkZGluZ3MgZm9yIGVhY2ggc3RhdGUKCmBgYHtyfQptdV9oYXQgPSBwaGlfeFssIGxhcHBseSguU0QsIG1lYW4pLCAuU0Rjb2xzID0gbmFtZXMocGhpX3gpWzI6bmNvbChwaGlfeCldLCBieSA9IHN0YXRlXVtvcmRlcihzdGF0ZSldCm11X2hhdApgYGAKCiMjIyBNb2RlbCAxIC0gRnJlcXVlbnRpc3QKTm9ybWFsIExBU1NPIHJlZ3Jlc3Npb24gb24gdGhlIGtlcm5lbCBtZWFuIGVtYmVkZGluZ3MgLSBubyBCYXllc2lhbiB0cmVhdG1lbnQuCgokJHlfaSA9IFxiZXRhXlRcaGF0e1xtdX1faSArIGIkJApEb2Vzbid0IGFjY291bnQgZm9yIHVuY2VydGFpbnR5IGluIHRoZSBjb2VmZmljaWVudHMgJFxiZXRhJCBvciBpbiB0aGUgbWVhbiBlbWJlZGRpbmdzICRcaGF0e1xtdX1faSQsIG9yIG9ic2VydmF0aW9uIGVycm9yCgojIyMjIEJ1aWxkIHRoZSBtb2RlbAoKVXNlIDUtZm9sZCBDViBrZXJuZWwgTEFTU08gdG8gZXN0aW1hdGUgJFxiZXRhJCwgJGIkIGZvciAkXGhhdHt5fV9pID0gXGJldGFeVFxoYXR7XG11fV9pICsgYiQuICBVc2UgdGhlIGZpcnN0IExBU1NPIHRvIHNlbGVjdCBub24temVybyBjb2VmcywgdGhlbiByZS1maXQgd2l0aG91dCBwZW5hbHR5IHRvIGFjY29udCBmb3IgY292YXJpYXRlIHNocmlua2FnZS4KCmBgYHtyfQojIGZpdCBpbml0aWFsIGxhc3NvIHRvIGZpbmQgc2lnbmlmaWNhbnQgY292YXJzCmZpdF9sYW1iZGEgPSBjdi5nbG1uZXQoeCA9IGFzLm1hdHJpeChtdV9oYXRbLCAtMSwgd2l0aCA9IEZdKQogICAgICAgICAgICAgICAgICAgICAgICwgeSA9IFlfdHJhaW5fYWdnW29yZGVyKHN0YXRlKSwgXSR5X2RlbV9wY3QKICAgICAgICAgICAgICAgICAgICAgICAsIG5mb2xkcyA9IDcpCgojIGdldCBpbmRpY2llcyBvZiBub24tMCBjb3ZhcnMgKGV4Y2VwdCBpbnRlcmNlcHQpCmluZCA9IHdoaWNoKGNvZWYoZml0X2xhbWJkYSwgcyA9ICdsYW1iZGEubWluJylbLTFdICE9IDApCgojIHJlLWZpdCBMQVNTTyBvbmx5IHdpdGggbm9uLTAgY292YXJzCmZpdCA9IGdsbW5ldCh4ID0gYXMubWF0cml4KG11X2hhdFssIGluZCArIDEsIHdpdGggPSBGXSkKICAgICAgICAgICAgICwgeSA9IFlfdHJhaW5fYWdnW29yZGVyKHN0YXRlKSwgXSR5X2RlbV9wY3QpCmBgYAoKCiMjIyMgUHJlZGljdAoKR2V0IGVtYmVkZGluZ3Mgb2YgdGVzdCBwb2ludHMgYW5kIGVtcGlyaWNhbCBtZWFucyBvZiBiYWdzIG9mIGVtYmVkZGVkIHRlc3QgcG9pbnRzLiAgVXNlIHNhbWUgbGFuZG1hcmsgcG9pbnRzICR1JCBhcyBmb3IgdGhlIHRyYWluaW5nIHNldC4KYGBge3J9CiMgZ2V0IGVtYmVkZGluZ3Mgb2YgdGVzdCBwb2ludHMKcGhpX3hfdGVzdCA9IGtlcm5lbE1hdHJpeChyYmYsIHggPSBhcy5tYXRyaXgoWF90ZXN0WywgLTFdKSwgeSA9IHUpCnBoaV94X3Rlc3QgPSBkYXRhLnRhYmxlKHN0YXRlID0gWF90ZXN0JHN0YXRlLCBwaGlfeF90ZXN0KQpzZXRuYW1lcyhwaGlfeF90ZXN0LCBjKCdzdGF0ZScsIHBhc3RlMCgndScsIDE6bl9jZW50ZXJzKSkpCgojIGNhbGN1bGF0ZSBlbXBpcmljYWwgbWVhbnMKbXVfaGF0X3Rlc3QgPSBwaGlfeF90ZXN0WywgbGFwcGx5KC5TRCwgbWVhbiksIC5TRGNvbHMgPSBuYW1lcyhwaGlfeF90ZXN0KVsyOm5jb2wocGhpX3hfdGVzdCldLCBieSA9IHN0YXRlXVtvcmRlcihzdGF0ZSldCm11X2hhdF90ZXN0X2dlbmRlciA9IGNiaW5kKHBoaV94X3Rlc3QsIHNleF9mZW1hbGUgPSBhcy5udW1lcmljKFhfdGVzdCRzZXhfZmVtYWxlID4gMCkpWywgbGFwcGx5KC5TRCwgbWVhbiksIC5TRGNvbHMgPSBuYW1lcyhwaGlfeF90ZXN0KVsyOm5jb2wocGhpX3hfdGVzdCldLCBieSA9IC4oc3RhdGUsIHNleF9mZW1hbGUpXVtvcmRlcihzdGF0ZSwgc2V4X2ZlbWFsZSldCmBgYAoKCkZpcnN0LCB0ZXN0IHByZWRpY3Rpb25zIGF0IHRoZSBzdGF0ZS1sZXZlbC4gIFdlIGNhbiBhbHNvIHRlc3QgYXQgdGhlIHN0YXRlLWdlbmRlciBsZXZlbC4KYGBge3J9CnlfaGF0X2FnZyA9IHByZWRpY3QoZml0LCBuZXd4ID0gYXMubWF0cml4KG11X2hhdF90ZXN0WywgaW5kICsgMSwgd2l0aCA9IEZdKSwgcyA9IDApCnlfaGF0X2FnZyA9IGNiaW5kKFlfdGVzdFssIC4oeV9kZW1fcGN0ID0gbWVhbih5X2RlbSksIC5OKSwgYnkgPSBzdGF0ZV1bb3JkZXIoc3RhdGUpXSwgeV9oYXRfYWdnKQoKZ2dwbG90KHlfaGF0X2FnZywgYWVzKHggPSB5X2RlbV9wY3QsIHkgPSBgMWAsIHNpemUgPSBOKSkgKyAKICBnZW9tX3BvaW50KCkgKyAKICB5bGFiKCJQcmVkaWN0ZWQgJSBEZW0iKSArCiAgeGxhYigiQWN0dWFsICUgRGVtIikgKwogIGdndGl0bGUoIlByZWRpY3RlZCBieSBhY3R1YWwgJSBEZW0gc3VwcG9ydCBhY3Jvc3Mgc3RhdGVzIikKCnlfaGF0X2FnZ19nZW5kZXIgPSBwcmVkaWN0KGZpdCwgbmV3eCA9IGFzLm1hdHJpeChtdV9oYXRfdGVzdF9nZW5kZXJbLCBpbmQgKyAxLCB3aXRoID0gRl0pLCBzID0gMCkKeV9oYXRfYWdnX2dlbmRlciA9IGNiaW5kKGNiaW5kKFlfdGVzdCwgc2V4X2ZlbWFsZSA9IGFzLm51bWVyaWMoWF90ZXN0JHNleF9mZW1hbGUgPiAwKSlbLCAuKHlfZGVtX3BjdCA9IG1lYW4oeV9kZW0pLCAuTiksIGJ5ID0gLihzdGF0ZSwgc2V4X2ZlbWFsZSldW29yZGVyKHN0YXRlLCBzZXhfZmVtYWxlKV0sIHlfaGF0X2FnZ19nZW5kZXIpCgpnZ3Bsb3QoeV9oYXRfYWdnX2dlbmRlciwgYWVzKHggPSB5X2RlbV9wY3QsIHkgPSBgMWAsIHNpemUgPSBOLCBjb2xvciA9IGFzLmZhY3RvcihzZXhfZmVtYWxlKSkpICsgCiAgZ2VvbV9wb2ludCgpICsgCiAgZ2VvbV9hYmxpbmUoaW50ZXJjZXB0ID0gMCwgc2xvcGUgPSAxKSArCiAgeWxhYigiUHJlZGljdGVkICUgRGVtIikgKwogIHhsYWIoIkFjdHVhbCAlIERlbSIpICsKICBnZ3RpdGxlKCJQcmVkaWN0ZWQgYnkgYWN0dWFsICUgRGVtIHN1cHBvcnQgYWNyb3NzIHN0YXRlcyAmIGdlbmRlcnMiKQoKcHJpbnQocGFzdGUwKCJNU0UgKHN0YXRlKTogIiwgbWVhbigoeV9oYXRfYWdnWyxgMWBdIC0geV9oYXRfYWdnJHlfZGVtX3BjdCleMikpKQpwcmludChwYXN0ZTAoIk1TRSAoc3RhdGUgJiBnZW5kZXIpOiAiLCBtZWFuKCh5X2hhdF9hZ2dfZ2VuZGVyWyxgMWBdIC0geV9oYXRfYWdnX2dlbmRlciR5X2RlbV9wY3QpXjIpKSkKCmBgYAoKUHJlZGljdCBhdCB0aGUgaW5kaXZpZHVhbC1sZXZlbApgYGB7cn0KIyBwcmVkaWN0IG9uIG5ldyAoaW5kaXZpZHVhbC1sZXZlbCkgZGF0YSBzZXQgb2YgWF90ZXN0IGVtYmVkZGVkIGluIGZlYXR1cmUgc3BhY2UKeV9oYXQgPSBwcmVkaWN0KGZpdCwgbmV3eCA9IGFzLm1hdHJpeChwaGlfeF90ZXN0WywgaW5kICsgMSwgd2l0aCA9IEZdKSwgcz0gMCkgICNubyByZWd1bGFyaXphdGlvbgp5X2hhdCA9IGNiaW5kKFlfdGVzdCwgeV9oYXQpCgojIGNhbGN1bGF0ZSBvdmVyYWxsIGNsYXNzaWZpY2F0aW9uIHJhdGUKY2xhc3NfcmF0ZXMgPSB5X2hhdFssIC4oLk4KICAgICAgICAgICwgcGN0X2NvcnJlY3QgPSBtZWFuKChgMWAgPiAwLjUgJiB5X2RlbSA9PSAxKSB8IChgMWAgPD0gMC41ICYgeV9kZW0gPT0gMCkpCiAgICAgICAgICAsIHBjdF9mYWxzZV9uZWcgPSBtZWFuKGAxYCA8PSAwLjUgJiB5X2RlbSAhPSAwKQogICAgICAgICAgLCBwY3RfZmFsc2VfcG9zID0gbWVhbihgMWAgPiAwLjUgJiB5X2RlbSAhPSAxKQogICAgICAgICAgLCBtc2UgPSBtZWFuKCh5X2RlbSAtIGAxYCleMikKICAgICAgICAgICldCmNsYXNzX3JhdGVzCmBgYAoKUGxvdCBhY3R1YWwgcGVyY2VudCBEZW0gYnkgZGVjaWxlIG9mIHRoZSBzY29yZQpgYGB7cn0KIyBnZXQgZGVjaWxlcyBvZiBwcmVkaWN0ZWQgc2NvcmUKeV9oYXRbLCBwcmVkX2RlY2lsZSA6PSBjdXQoYDFgLCBicmVha3MgPSBxdWFudGlsZShgMWAsIHByb2JzID0gc2VxKDAsMSwwLjEpKSwgbGFiZWxzID0gMToxMCwgaW5jbHVkZS5sb3dlc3QgPSBUKV0KCiNnZ3Bsb3QoeV9oYXQsIGFlcyh4ID0gYDFgLCBjb2xvciA9IGFzLmZhY3Rvcih5X2RlbSkpLCBhbHBoYSA9IDAuMikgKyBnZW9tX2RlbnNpdHkoKQoKI3Bsb3QgYWN0dWFsIHBjdCBvZiBlYWNoIGRlY2lsZSB0aGF0IHN1cHBvcnRzIGRlbXMKZ2dwbG90KHlfaGF0WywgLihwY3RfeV9kZW0gPSBtZWFuKHlfZGVtKSksIGJ5ID0gcHJlZF9kZWNpbGVdLCBhZXMoeCA9IHByZWRfZGVjaWxlLCB5ID0gcGN0X3lfZGVtKSkgKyAKICBnZW9tX2JhcihzdGF0ID0gJ2lkZW50aXR5JykgKwogIGdndGl0bGUoIlBjdCBEZW0gc3VwcG9ydGVyIGJ5IHNjb3JlIGRlY2lsZSIpCmBgYAoKCgojIyBNb2RlbCAyIC0gQmF5ZXNpYW4gTGluZWFyIFJlZ3Jlc3Npb24KCkFkZGluZyB1bmNlcnRhaW50eTogaW5jb3Jwb3JhdGUgdW5jZXJ0YWludHkgb2YgdGhlIGNvZWZmaWNpZW50cyAkXGJldGEkCgpXZSBjYW4gaW50ZXJwcmV0IHRoaXMgaW4gMiB3YXlzOiAxKSBCYXllc2lhbiBsaW5lYXIgcmVncmVzc2lvbiAyKSBHYXVzc2lhbiBwcm9jZXNzIHJlZ3Jlc3Npb24gY29uZGl0aW9uZWQgb24gZmluaXRlIHNldCBvZiBvYnNlcnZlZCBwb2ludHMuCgoxKSBCYXllc2lhbiBsaW5lYXIgcmVncmVzc2lvbgokJFxiZXRhIFxzaW0gXG1hdGhjYWx7Tn0oMCwgXFNpZ21hX3ApXFwKeV9pIHwgXGJldGEsIFxtYXRoYmZ7eH1faSBcc2ltIFxtYXRoY2Fse059KFxiZXRhXlRcaGF0e1xtdX1faSwgXHNpZ21hXjIpJCQKCkZvciBkaXNpdHJpYnV0aW9uIHJlZ3Jlc3Npb24sIHdlIHN1YnN0aXR1dGUgaW4gdGhlIGVtcGlyaWNhbCBtZWFuIGVtYmVkZGluZ3MgJFxoYXR7XG11fV9pJCBmb3IgJHhfaSQuCgoyKSBHUCByZWdyZXNzaW9uIGNvbmRpdGlvbmVkIG9uIG9ic2VydmVkIHBvaW50cyAod2VpZ2h0IHNwYWNlIHZpZXcpCiQkCmYoXG1hdGhiZnt4fSkgPSBccGhpKFxtYXRoYmZ7eH0pXlRcbWF0aGJme3d9XFwKeSA9IGYoXG1hdGhiZnt4fSkgKyBcZXBzaWxvblxcClxlcHNpbG9uIFxzaW0gXG1hdGhjYWx7Tn0oMCwgXHNpZ21hX25eMilcXApcbWF0aGJme3d9IFxzaW0gXG1hdGhjYWx7Tn0oMCwgXFNpZ21hX3ApCiQkCgp3ZSBjYW4gZXF1aXZhbGVudGx5IHdyaXRlCiQkCmYgXHNpbSBcbWF0aGNhbHtHUH0oMCwgayh4LCB4JykpXFwKeV9pfGYoeF9pKSBcc2ltIFxtYXRoY2Fse059KDAsIEsgKyBJXHNpZ21hXjIpCiQkCndoZXJlICRrKHgseCcpJCBpcyB0aGUga2VybmVsIGNvcnJlc3BvbmRpbmcgdG8gdGhlIGZlYXR1cmUgbWFwICRccGhpJAoKRm9yIEdQUiwgdGhlIGZlYXR1cmUgbWFwICRccGhpKHhfaSkkIGNvcnJlc3BvbmRzIHRvIHRoZSBrZXJuZWwgdGhhdCB3YXMgdXNlZCB0byBlbWJlZCAkeF9pJCBpbiBmZWF0dXJlIHNwYWNlLiAgVXNpbmcgbGFuZG1hcmsgcG9pbnRzICRcbWF0aGJme3V9ID0gXHt1X2xcfV97bD0xfV5kJCwgdGhpcyBiZWNvbWVzICRccGhpKHhfaSkgPSBbayh4X2ksIHVfMSksIFxkb3RzLCBrKHhfaSwgdV9sKV0kCgpXZSBjYW4gaW50ZXJwcmV0IHRoZSB3ZWlnaHRzICRcbWF0aGJme3d9JCBvZiBHUCByZWdyZXNzaW9uIGFzIHRoZSBjb2VmZmljaWVudHMgJFxiZXRhJCBpbiBCTFIuICAkXFNpZ21hX3AkIGlzIHRoZSBwcmlvciBjb3ZhcmlhbmNlIG1hdHJpeCBmb3IgdGhlIHdlaWdodHMvY29lZmZpY2llbnRzLgoKRml0IEdQCmBgYHtyfQojIG9ic2VydmF0aW9ucyA9IG11X2hhdAojIG91dGNvbWVzID0gWV90cmFpbl9hZ2ckeV9kZW1fcGN0CgpmaXRfZ3AgPSBnYXVzc3ByKHggPSBhcy5tYXRyaXgobXVfaGF0WywgLTFdKSwgeSA9IFlfdHJhaW5fYWdnW29yZGVyKHN0YXRlKSxdJHlfZGVtX3BjdCwga2VybmVsID0gJ3ZhbmlsbGFkb3QnKQpmaXRfZ3AKCiNmaXRfZ3BAYWxwaGEKYGBgCgpOb3csIHdlIHByZWRpY3QgdGhlIEdQIGF0IHRoZSBiYWctbGV2ZWwKYGBge3J9Cgp5X2hhdF9hZ2cgPSBwcmVkaWN0KGZpdCwgbmV3eCA9IGFzLm1hdHJpeChtdV9oYXRfdGVzdFssIGluZCArIDEsIHdpdGggPSBGXSksIHMgPSAwKQp5X2hhdF9hZ2cgPSBjYmluZChZX3Rlc3RbLCAuKHlfZGVtX3BjdCA9IG1lYW4oeV9kZW0pLCAuTiksIGJ5ID0gc3RhdGVdW29yZGVyKHN0YXRlKV0sIHlfaGF0X2FnZykKCmdncGxvdCh5X2hhdF9hZ2csIGFlcyh4ID0geV9kZW1fcGN0LCB5ID0gYDFgLCBzaXplID0gTikpICsgCiAgZ2VvbV9wb2ludCgpICsgCiAgeWxhYigiUHJlZGljdGVkICUgRGVtIikgKwogIHhsYWIoIkFjdHVhbCAlIERlbSIpICsKICBnZ3RpdGxlKCJQcmVkaWN0ZWQgYnkgYWN0dWFsICUgRGVtIHN1cHBvcnQgYWNyb3NzIHN0YXRlcyIpCmBgYAoKQW5kIGF0IHRoZSBpbmRpdmlkdWFsIGxldmVsCmBgYHtyfQp5X2hhdF9ncCA9IHByZWRpY3QoZml0X2dwLCBuZXdkYXRhID0gcGhpX3hfdGVzdFssIC0xXSkKeV9oYXRfZ3AgPSBjYmluZChZX3Rlc3QsIHlfaGF0X2dwKQoKeV9oYXRfZ3BbLCBwcmVkX2RlY2lsZSA6PSBjdXQoVjEsIGJyZWFrcyA9IHF1YW50aWxlKFYxLCBwcm9icyA9IHNlcSgwLDEsMC4xKSksIGxhYmVscyA9IDE6MTAsIGluY2x1ZGUubG93ZXN0ID0gVCldCgojcGxvdCBhY3R1YWwgcGN0IG9mIGVhY2ggZGVjaWxlIHRoYXQgc3VwcG9ydHMgZGVtcwpnZ3Bsb3QoeV9oYXRfZ3BbLCAuKHBjdF95X2RlbSA9IG1lYW4oeV9kZW0pKSwgYnkgPSBwcmVkX2RlY2lsZV0sIGFlcyh4ID0gcHJlZF9kZWNpbGUsIHkgPSBwY3RfeV9kZW0pKSArIAogIGdlb21fYmFyKHN0YXQgPSAnaWRlbnRpdHknKSArCiAgZ2d0aXRsZSgiUGN0IERlbSBzdXBwb3J0ZXIgYnkgc2NvcmUgZGVjaWxlIikKYGBgCgoKCgoKCgoKCiMjIFRISU5HUyBUTyBETwoKIlVubGFiZWxlZCBkYXRhIGNhbiBzaW1wbHkgYmUgdXNlZCBhcyBjZW50ZXJzIiAtIFN6YWJvIDIwMTYKLSBub3QgZW50aXJlbHkgYXBwbGljYWJsZSBoZXJlLCB3ZSBvYnNlcnZlIGxhYmVscywganVzdCBhdCB0aGUgYWdncmVnYXRlIGxldmVsCgoqR2VuZXJhbCBxdWVzdGlvbnMqCgotIGhvdyBtdWNoIHRvIGZvY3VzIG9uIHRoaXMgc3BlY2lmaWMgZXhhbXBsZSB2LiB0aGUgbW9yZSBnZW5lcmFsIHF1ZXN0aW9uIG9mICJkb2VzIGRpc3QgcmVnIHdvcmsgZm9yIGluZGl2aWR1YWxzIj8KLSBXaGF0IGFyZSB0aGUgcmVhc29ucyB0aGF0IHlvdSdyZSBza2VwdGljYWwgdGhhdCBCRFIgd2lsbCB3b3JrIGZvciBpbmRpdmlkdWFscz8gIFdoYXQgYXJlIHRoZSBzaG9ydGNvbWluZ3MgdGhhdCB5b3UgZW52aXNpb24/CgotIGFkdmFudGFnZXMgdG8gaW1wbGVtZW50aW5nIGJhc2UgbW9kZWwgYXMgR1BSIHYuIEJMUj8gIEdQUiBtb3JlIGZsZXhpYmxlIGZvciBub24tR2F1c3NpYW4gb3V0Y29tZXMgKGkuZS4gbXVsdGlub21pYWwpPyAgQ2FuIGNvbWJpbmUga2VybmVscyAoaS5lLiBtZWFuIGVtYmVkZGluZ3MgKyBzcGF0aWFsIGNvbXBvbmVudCk/ICBEb3duc2lkZSBpcyBjb21wdXRhdGlvbj8KCi0gYmF0Y2hpbmc/CgoqQ2hvb3NlIGRhdGEgc2V0KjoKCi0gY2Vuc3VzIG1pY3JvZGF0YSB3aXRoIGFjdHVhbCBjb3VudHktbGV2ZWwgZWxlY3Rpb24gb3V0Y29tZXMuIHBybzogbGFyZ2UgKGxpa2UgYWN0dWFsIHZvdGVyZmlsZSksIGNvbjogb3V0Y29tZXMgbm90IGF0IHRoZSBpbmRpdmlkdWFsLWxldmVsCi0gUGV3IHN1cnZleSBkYXRhLiBwcm86IG91dGNvbWVzIGF0IHRoZSBpbmRpdmlkdWFsLWxldmVsLCBjb246IGl0J3Mgc21hbGwKLSBjZW5zdXMgbWljcm9kYXRhIGZvciBjb3ZhciBkYXRhIHdpdGggc3VydmV5IGRhdGEgZm9yIG91dG9tZXMuICBQcm86IGxhcmdlIGNvdmFyIGRhdGEsIGNsb3NlciB0byByZWFsIHNjZW5hcmlvLCBjb246IGNhbiBvbmx5IGxpbmsgYXQgdGhlIHN0YXRlLWxldmVsLCBzdGlsbCBubyBpbmRpdmlkdWFsLWxldmVsIG91dGNvbWUgZGF0YQotIHN5bnRoZXRpYyBkYXRhLiBwcm86IGNhbiBtYWtlIGl0IGV4YWN0bHkgd2hhdCB3ZSB3YW50LCBjb246IGVmZm9ydCB0byBzeW50aGVzaXplCgoqKiBiZW5lZml0IG9mIHRoZSBsYXJnZXIgZGF0YSBzZXQgLSBtb3JlIHJlYWxpc3RpYywgd2UnZCByZWFsbHkgbmVlZCB0aGUgbGFuZG1hcmsgcG9pbnRzIHdoZW4gd2UgZG9uJ3QgaWYgb25seSB3b3JraW5nIHdpdGggdGhlIHN1cnZleSBkYXRhCgoqKiBsYXJnZXIgZGF0YSBhbHNvIGhhcyBtb3JlIGNvdmFyaWF0ZSBkYXRhIGF2YWlsYWJsZQoKCipwYXJhbWV0ZXIgdHVuaW5nKjoKCi0gbnVtYmVyIG9mIGxhbmRtYXJrIHBvaW50cyAoaG93IGlzIHRoaXMgbm9ybWFsbHkgY2hvc2VuPykgLSBMZW9uJ3MgY29kZSBoYXMgMzAgYXMgdGhlIGRlZmF1bHQgZm9yIDEwMDAgdGVzdCwgMTAwMCB0cmFpbmluZyAoc2ltaWxhciBzY2FsZSB0byBkYXRhIGhlcmUpCi0gQ2FsY3VsYXRlIGxhbmRtYXJrIHBvaW50cyB3aXRoIGZ1bGwgZGF0YSBzZXQgb3IganVzdCB0cmFpbmluZz8KLSBzY2FsZSBwYXJhbWV0ZXIgJFxzaWdtYSQgZm9yIFJCRiBrZXJuZWwgKHJpZ2h0IG5vdyBqdXN0IHVzaWluZyBtZWRpYW4gaGV1cmlzdGljKQoKKkNob29zaW5nIGJhZyBkZWZpbml0aW9ucyoKCi0gU2luY2Ugd2UgaGF2ZSBpbmRpdmlkdWFsLWxldmVsIGxhYmVscyBpbiB0aGUgc3VydmV5IGRhdGEsIHdlIGNhbiBjaG9vc2UgaG93IHdlIGFnZ3JlZ2F0ZSBpdCAoZG8gd2UgZXZlbiBuZWVkIHRvIGFnZ3JlZ2F0ZSBpdD8gWUVTIGJlY2F1c2Ugd2UgbmVlZCB0aGUgc3RyYXRhIGRlZmluaXRpb25zIGluIG9yZGVyIHRvIGxpbmsgaXQgYmFjayB0byB0aGUgZnVsbCBmaWxlKQotIEJhZ3MgYmFzZWQgb24gZGVtb3MgaW5zdGVhZCBvZiBzdGF0ZSAobW9yZSBzcGVjaWZpYyB0byB0aGlzIHBvbGwgdGhhbiBtb3N0IHdlJ2QgYmUgZGVhbGluZyB3aXRoKT8KCipSQkYgTmV0d29yayB2LiBTaW5nbGUgUkJGIGtlcm5lbCoKCi0gcmVsYXRpb25zaGlwIGJldHdlZW4gbmV1cmFsIG5ldHdvcmsgaW1wbGVtZW50YXRpb24gKHJlIExlb24ncyBjb2RlKSBhbmQgdGhlIGtlcm5lbCBtZWFuIGVtYmVkZGluZz8gSnVzdCB0aGF0IHdlIGNhbiBleGFjdGx5IGludGVycHJldCBrZXJuZWwgZW1iZWRkaW5nIGFzIGEgbmV1cmFsIG5ldHdvcmsgYW5kIG1heWJlIGVhc2llciB0byBpbXBsZW1lbnQvZmFzdCBvcHRpbWl6YXRpb24gb2YgbWFyZ2luYWwgbG9nIGxpa2VsaWhvb2Q/CgoqQmVuY2htYXJrKjogCgotIENvbXBhcmUgYWdhaW5zdCBsYXNzby9jdXN0b20tYnVpbHQgbG9naXQgb3IgbWF5YmUgTk4gdGhhdCBoYXMgYWNjZXNzIHRvIHRoZSBpbmRpdmlkdWFsLWxldmVsIG91dGNvbWVzIChtb3JlIHJlYWxpc3RpYyBiZW5jaG1hcmsgdGhhbiksIGp1c3QgdXNpbmcgYSBiaWFzZWQgc3Vic2V0IG9mIHRoZSB0cmFpbmluZyBvYnNlcnZhdGlvbnMgKG1pbWljIGFzeW1tZXRyaWMgbWF0Y2hpbmcpCi0gTVJQIC0gc3RhbmRhcmQgZm9yIEVJIGluIHBvbGl0aWNhbCBhcHBsaWNhdGlvbnMKLSBzdHJhaWdodCBzdHJhdGlmaWNhdGlvbiB1c2luZyBjb3ZhcmlhdGVzIG9ic2VydmVkIGluIHRoZSBzdXJ2ZXkKCgpJbXBsZW1lbnQgKmZ1bGwgQmF5ZXNpYW4gdHJlYXRtZW50KiAodW5jZXJ0YWludHkgaW4gbWVhbiBlbWJlZGRpbmdzIGFzIHdlbGwgYXMgYmFnIHNpemUpCgotIGZsZXNoIG91dCBmdWxsIGdlbmVyYXRpdmUgbW9kZWwKCipTdHJ1Y3R1cmVkIG9ic2VydmF0aW9uYWwgdW5jZXJ0YWludHkqCgotIGhlcmlhcmNoaWNhbCBtb2RlbD8KCgoKCgoKCgoKCg==